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Gravitational waves from oscillating neutron stars in axial symmetry are 
studied performing numerical simulations in full general relativity. Neutron 
' stars are modeled by a polytropic equation of state for simplicity. A gauge- 
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invariant wave extraction method as well as a quadrupole formula are adopted 
for computation of gravitational waves. It is found that the gauge- invariant 
variables systematically contain numerical errors generated near the outer 
boundaries in the present axisymmetric computation. We clarify their origin, 
and illustrate it possible to eliminate the dominant part of the systematic 
errors. The best corrected waveforms for oscillating and rotating stars cur- 



^ ■ rently contain errors of magnitude ~ 10 ^ in the local wave zone. Comparing 

the waveforms obtained by the gauge-invariant technique with those by the 
quadrupole formula, it is shown that the quadrupole formula yields approx- 
imate gravitational waveforms besides a systematic underestimation of the 
amplitude of 0{M/R) where M and R denote the mass and the radius of 
neutron stars. However, the wave phase and modulation of the amplitude 
can be computed accurately. This indicates that the quadrupole formula is 
a useful tool for studying gravitational waves from rotating stellar core col- 
lapse to a neutron star in fully general relativistic simulations. Properties of 
the gravitational waveforms from the oscillating and rigidly rotating neutron 
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stars are also addressed paying attention to the oscillation associated with 

fundamental modes. 
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I. INTRODUCTION 



One of the most important roles of numerical simulations in general relativity is to 
predict gravitational waveforms emitted by general relativistic and dynamical astrophysical 
phenomena. Rotating stellar core collapse and nonspherical oscillation of neutron stars 
are among the possible sources of gravitational waves. Therefore, fully general relativistic 
numerical simulation for them is an important subject in this field [1] . 

To date, there has been no systematic work for computation of gravitational waves 
from rotating stellar core collapse to a neutron star in fully general relativistic simulation 
(but see [2]). The gravitational waveforms have been computed only in the Newtonian 
gravity [3-9] or in an approximate general relativistic gravity [1] using the so-called conformal 
flatness approximation (or Isenberg- Wilson- Mathews approximation). As demonstrated in 
[1], general relativistic effects modify the evolution of the collapse and emitted gravitational 
waveforms significantly. Thus, the simulation in full general relativity appears to be the 
best approach for accurate computation of gravitational waves. 

In the case that the progenitor of the core collapse is not very rapidly rotating, nonax- 
isymmetric instabilities do not set in and, hence, the collapse will proceed in an axisymmetric 
manner. In such a collapse, the amplitude of gravitational waves measured in a local wave 
zone at r ~ A where A denotes the gravitational wave length will be smaller, by two or 
three orders of magnitude, than that in highly nonaxisymmetric phenomena such as merg- 
ers of binary neutron stars and black holes. The amplitude of gravitational waves from an 
oscillating neutron star is also likely to be small due to its small nonspherical deformation. 
Technically, it is not easy to extract gravitational waves of small amplitude from metric 
computed in numerical simulations, in which a numerical noise is in general contained. The 
numerical noise is generated due to the following reasons: 

Gravitational waves are usually extracted from the metric in the wave zone in general 
relativistic simulations. Although they should be extracted at the null infinity, the outer 
boundaries of computational domain are located at a finite radius whenever the 3-1-1 for- 
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malisms are adopted. Thus, the outer boundary conditions are imposed at finite radii and 
in general they are approximate conditions. As a result, a small numerical error may be ex- 
cited around the outer boundaries. Here, the possible candidates of the numerical error are 
unphysical nonwave modes, spurious gauge modes, back reflections at the outer boundaries, 
and roundoff errors. 

In this paper, we study gravitational waves from oscillating neutron stars in axial sym- 
metry. Neutron stars in equilibrium are simply modeled by n = 1 polytropes. Oscillations 
of neutron stars are followed by axisymmetric numerical simulations in full general relativ- 
ity. Gravitational waves are extracted using a gauge-invariant wave extraction technique. 
The gauge-invariant variables are not contaminated by gauge modes and, hence, we can 
focus on other error sources using this variables. We also adopt a quadrupole formula for 
approximately computing gravitational waveforms to clarify its validity. 

This work was planned from the following four motivations. The flrst one is to specify 
the error sources contained in the gauge-invariant variables extracted in the local wave zone. 
As mentioned above, they could be contaminated by nonwave components and numerical 
errors. In particular, it is important to specify systematic error components contained in 
the gauge-invariant variables since as indicated in Sec. IV, the systematic errors may be 
eliminated at least partly if their origin is clarified. 

The second motivation is to understand how large computational domains are needed to 
extract gravitational waveforms within ~ 10% error. Since the gauge-invariant variables are 
extracted at finite radii, gravitational waveforms (in particular the amplitude) are slightly 
different from the asymptotic ones. It is important to clarify how magnitude of the error 
depends on the radius at which we impose the outer boundary conditions and on the radius 
at which we extract gravitational waves. A similar study was carried out about 15 years ago 
by Abrahams and Evans [10]. However, they were interested only in specific gauge conditions 
which were often used in axisymmetric numerical simulations in general relativity at that 
time. Moreover, the simulations were carried out only for non-rotating stars. In this paper, 
we adopt a different gauge condition often used nowadays in three-dimensional simulations, 
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and report numerical results both for nonrotating and rotating stars. 

The third motivation, in which we are most interested in the present study, is to investi- 
gate validity of a quadrupole formula in fully general relativistic simulations. For computa- 
tion of gravitational waves generated by oscillations of gravitational field such as quasinormal 
mode ringings of black holes, quadrupole formulas cannot work. However, in rotating stellar 
core collapse to a neutron star and in oscillating neutron stars in which gravitational waves 
are generated mainly by matter motions, quadrupole formulas may be able to yield an accu- 
rate waveform. This method can be applied much more easily than geometrical methods in 
which gravitational waves are extracted from metric in the wave zone. Thus, a quadrupole 
formula which can yield high-quality approximate waveforms will be a robust method for 
computing gravitational waves of small amplitude from a noisy numerical data set. Note 
that a similar work has been already done by Siebel et al. [11,2] in a null-cone formulation. 
We here carry out the similar study for a 3-1-1 approach. 

The last motivation is to understand the properties of oscillations of rotating neutron 
stars. During rotating stellar core collapse, gravitational waves associated with oscillations 
of a formed protoneutron star are hkely to be emitted (see, e.g., [1]). Prom the study for 
oscillating and rotating neutron stars, we will be able to understand what oscillation modes 
are relevant for the emission of gravitational waves during core collapses. We here pay 
attention to two fundamental oscillation modes (quasiradial and quadrupole p modes of no 
node for the density perturbation) which are candidates for the dominant modes in the 
oscillating and rotating stars formed after the collapse. 

This paper is organized as follows. In Sec. II, our numerical implementations for axisym- 
metric general relativistic simulation are briefly reviewed. In Sec. Ill, the gauge-invariant 
wave extraction technique and the quadrupole formula adopted in the present work are 
described. Sec. IV presents numerical results of gravitational waveforms emitted from os- 
cillating neutron stars. The simulations were performed both for nonrotating and rotating 
neutron stars using an axisymmetric code recently developed [12]. Sec. IV is devoted to a 
summary. Throughout this paper, we adopt the geometrical units in which G — c — 1 where 
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G and c are the gravitational constant and the speed of hght, respectively. 



II. NUMERICAL IMPLEMENTATION 

A. Summary of formulation 

We performed fully general relativistic simulations in axial symmetry using the same 
formulation as that in [12], to which the readers may refer for details and basic equations. 
The fundamental variables for hydrodynamics are: 

p: rest mass density, 
e: specific internal energy, 
P: pressure, 
u'^: four velocity, 

where subscripts i, j, /c, • • • denote x, y, and z, and ^ the spacetime components. In addition, 
we define a weighted density p*(= pau^e^'^) and a weighted four- velocity (1+e+P/ p)ui). 
From these variables, the total baryon rest-mass and angular momentum of system, which 
are conserved quantities in axisymmetric spacetimes, can be defined as 



t*^ j d^xp^, (2) 
J = j d?xp^u^. (3) 



General relativistic hydrodynamic equations are solved using a so-called high-resolution 
shock-capturing scheme [14,15,12] in cylindrical coordinates (or on y = plane in Cartesian 
coordinates) . 

The fundamental variables for geometry are: 

a: lapse function, 
(3^: shift vector, 
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7ij: metric in 3D spatial hypersurface, 



7= e 



det(7ij), 



i^^i,: extrinsic curvature. 



(4) 



We evolve 7jj, 0, = e "''^(-ft^ij — 'jijKf^^), and the trace of the extrinsic curvature Kj} 
together with three auxiliary functions Fj = S^'^djjik with an unconstrained free evolution 
code as done in [16,17,19,20,12]. 

The Einstein equations are solved in the Cartesian coordinates. To impose axisymmetric 
boundary conditions, the Cartoon method is adopted [13]: Assuming a reflection symmetry 
with respect to the equatorial plane, we perform simulations using a fixed uniform grid with 
the grid size N x 3 x N for (x, y, z) which covers a computational domain as < a; < L, 
< ^ < L, and —A < y < A. Here, N and L are constants and A — L/N . For y — ±A, 
the axisymmetric boundary conditions are imposed using data sets on the y = plane. 

The slicing and spatial gauge conditions adopted in this paper are basically the same 
as those in [16-19]; i.e., we impose an "approximately" maximal slice condition {Kj^ ^ 0) 
and an "approximately" minimal distortion gauge condition [Di{dt'y''^) ~ where Di is the 
covariant derivative with respect to ^ij] [16,17,19]. In the approximately minimal distortion 
gauge condition, Fj is zero in the linear order in hij{= jij — Sij) if Fj = initially. Thus, in 
the wave zone, hij approximately satisfies a transverse condition hijj — 0. 

We also performed a few simulations using a dynamical gauge condition [21] in which we 
solve 



where At denotes a timestep in a numerical computation. We have found that the magnitude 
of the numerical error depends weakly on the spatial gauge condition, but gravitational 
waveforms and qualitative nature of the numerical error do not. Thus, in this paper, we 
present the results obtained in the approximately minimal distortion gauge condition. 



kf3' = ^'\Fi + AtdtFi), 



(5) 
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During numerical simulations, violations of the Hamiltonian constraint and conservation 
of mass and angular momentum are monitored as code checks. Numerical results for several 
test calculations, including stability and collapse of nonrotating and rotating neutron stars, 
have been described in [12]. Several convergence tests have been also presented in [12]. 



B. Outer boundary conditions 

Outer boundary conditions for geometric variables have been modified from previous 
ones [16-19]. We impose a boundary condition for cj) as 

0=^ + O(r-^), (6) 

where M denotes the ADM mass of a system computed a,t t — 0. We note that M is an 
approximately conserved quantity in axial symmetry, since only a small amount of gravita- 
tional waves are emitted in axisymmetric oscillations. 

For hij{— jij — Sij), we first carry out a coordinate transformation from the Cartesian 
coordinates to spherical polar coordinates (r, 9, (p) with the flat metric, rjab [the subscripts 
a and b denote components of the spherical polar coordinates {r,9,ip)], and then impose 
outgoing- wave outer boundary conditions of the form 

hffT^ = fi{t-r), 
heer^Mt-r), 

K^r ^ feit-r), (7) 

where /i-g denotes a tetrad component and /j (i = 1 ~ 6) denote functions: Since fi{t — r) — 
fi[t — At — (r — At)], the value of fi{t — r) at the outer boundaries should be equal to 
that at a time t — At and at a radius r — At which is determined using the values of 8 
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nearby grid points at a previous timestep [20] . These sets of the outer boundary conditions 
are well-suited for a solution of the linearized Einstein equation in the transverse-traceless 
gauge condition for hab if they are imposed in a distant wave zone [24] . In the local wave 
zone, however, equations (7) are approximate boundary conditions since higher order terms 
of 1/r are neglected. Therefore, numerical errors and unphysical solutions may be generated 
around the outer boundaries. In addition, hij should physically contain nonwave modes 
(such as stationary multipole modes) which do not obey the boundary conditions (7). Since 
no attention is paid to such modes in imposing the condition (7), additional numerical errors 
may be excited. 

In the case of the approximately minimal distortion gauge, we adopt the same boundary 
conditions for Fj, K, and A^j as those used in previous papers [16,17]: i.e., Fi — K — 0, and 
an outgoing- wave boundary condition is imposed for A^j. In the dynamical gauge condition, 
an outgoing- wave boundary condition is imposed for Fj, because it obeys a hyperbolic- type 
equation. 



Gravitational waves are extracted in terms of a gauge- invariant technique [22,23]. In this 
method, the fully nonlinear 3-metric 'jab in spherical polar coordinates is split as rjab + ^ab, 
where ^ab is regarded as the perturbation on the flat background. In axial symmetry, ^ab 
can be decomposed as J2i CibJ where Cib is given by 



III. WAVE EXTRACTION METHODS 



A. Gauge-invariant technique 



H21Y10 hiiYiofi 







\ 



* 



r\KiYio + GiWio) 




r"" sin^ 9{KiYio - GiWio) j 
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'^O CideYiosmO ^ 



+ 



* -r^DiWiQsmO 



* * 







(8) 



Here, * denotes the symmetric components. The first term in Eq. (8) corresponds to even 
parity (polar) perturbations and the second one to odd parity (axial) perturbations. The 
quantities H211 hu, Ki, Gi, Ci, and Di are functions of r and t, and are calculated by 
performing integrations over a two-sphere of a given radius. Yiq is the spherical harmonic 
function, and Wiq is defined as 



W, 
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{def-coiOde 



Yio. 



(9) 



\ 



The gauge-invariant variables of even and odd parities are then defined as [22,23] 

^^^^'"^^^ ^]^{2^2; + ^(^ + 1)^h}. (10) 

(11) 



where 



kn= Ki + l{l + l)Gi + 2rdrGi - 2^, 

r 

k^r= H„ - S- \r{Ki + /(/ + l)Gi} . 



dr 



(12) 
(13) 



Luminosity of gravitational waves is computed from 



dE 



dt 32n 



j:[\dtRfn\dtR?\'\ 



(14) 



The time derivative of the gauge- invariant quantities in Eq. (14) is taken by a finite- 
differencing. Hereafter, we focus only on the even-parity mode with / = 2 because for 
the oscillations considered in this paper, its amplitude is much larger than that of other 
modes. 

In an appropriate radiation gauge, -|— mode of gravitational waves is extracted from the 
asymptotic behavior of the following quantity at r — > 00; 
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For 1 — 2, can be written as 

where A2{t) denotes a function of time 
as 

pE 

B. Quadrupole formula 

In quadrupole formulas, gravitational waves are computed from 

K=P.^P,f,^), (18) 

where and Pi'' — 5ij — niUj {rii — x^r) denote a tracefree quadrupole moment and a 
projection tensor. Prom this expression, +-mode of gravitational waves with Z = 2 in axial 
symmetry is written as 

r 

where lij denotes an appropriately defined quadrupole moment, and lij its second time 
derivative. Equation (19) implies that in quadrupole formulas, A2{t) — Ixxitret) ~ hzitret)- 
tret is a retarded time which is approximately defined by 

tret = t - rcirc - 2M In - l) , (20) 

where Vdrc — r{l + M/2r)^ ^ M. Equation (20) is the vahd expression only for the 
Schwarzschild spacetime. In axisymmetric spacetimes, the retarded time should depend on 
the direction of wave propagation. However, the magnitude of the difference between Eq. 



^sin^^, (16) 
r 

Using A2, the asymptotic behavior of i?f is written 
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(20) and the exact one would be of 0{M) and, hence, much smaller than the typical wave 
length. 

In fully general relativistic and dynamical spacetimes, there is no unique definition for 
the quadrupole moment and, hence, for lij. Here, we choose for simplicity 

lij = j p^x'xU^x. (21) 

Then, using the continuity equation of the form 

9tp* + 9,(py)=0, (22) 

we can compute the first time derivative as 

hj - j p*{v'x^ + xV)d^x. (23) 

To compute I^j, we carried out a finite differencing of the numerical result for lij. 

Since the quadrupole moment 1^ is not defined uniquely in the dynamical spacetime in 
general relativity, gravitational waveforms computed by quadrupole formulas depend on the 
form of lij. In addition, they depend on the gauge conditions, since a physical point which 
coordinates x^ and t denote (and, as a result, the magnitude of lij and lij) is not identical in 
two different gauge conditions. Even if an identical definition of the quadrupole formula is 
employed, waveforms do not in general agree when different gauge conditions are adopted. 
Therefore, we should keep in mind that the waveforms computed from Eq. (19) are special 
ones obtained in specific choices of lij and the gauge condition. All these facts imply that 
to know the validity of the quadrupole formula which one chooses, comparison between the 
waveforms by the quadrupole formula with those extracted from metric should be made. 

IV. NUMERICAL RESULTS 
A. Setting 

The simulations were carried out along the following procedure: (1) neutron stars in 
equilibrium were prepared, (2) a perturbation to the equilibria was added, (3) the constraint 
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conditions were re-imposed by solving the constraint equations for the perturbed state, and 
then (4) we started the time evolution. 

To model neutron stars in equilibrium, we simply adopt the polytropic equation of state 

as 

P^Kp^+^, (24) 

where K is the polytropic constant and n the polytropic index. For the evolution of the 
neutron stars, we adopt a F-law equation of state in the form 

p = (r-i)p£, (25) 

where F = 1 + We set n = 1 (F = 2) as a reasonable qualitative approximation to cold, 
nuclear equation of state. With these equations of state, realistic neutron stars may not be 
modeled precisely. However, in the current situation that no one know a real equation of 
state of neutron stars exactly, modeling neutron stars with a simple equation of state is an 
adequate and popular strategy. 

In the present choice of the equation of state, physical units enter the problem only 
through the polytropic constant K, which can be chosen arbitrarily or else completely scaled 
out of the problem. For n — 1, K^^'^ has units of length, time, and mass, and units 
of density in the geometrical units. Using this property, we rescale all the quantities to be 
nondimensional and show only the nondimensional quantities. 

One often prefers to use particular dimensional units even in the polytropic equation of 
state. For example, in [15], the authors fix the value of K as 1.455 x 10^ cgs. For the sake of 
comparison with the previous paper, we convert non-dimensional quantities to dimensional 
ones with the polytropic constant K — 1.455 x 10^ cgs. In this special value, the mass, the 

density, and the time in the dimensional units are written as 

/ K \i/2/ M \ 



1.455 X 105 cgs/ V0.180 

Pdim = 1-86 X 10^^ g/cm^f ^ (—^], (27) 

^"^^ ^' V 1.455 X 105 cgs/ V0.300y' ^ ^ 

Tdi^ = 4.93msecf — • (28) 



A1.455 X 105 cgsy Viooy' 
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Pc 




M 


M/R 


J/M^ 


, //^ -T-iQ /O It r 1 /0\ 

Posc/(27ri?3/2M-l/2) 


SI 


0.127 


0.150 


0.140 


0.146 





0.84 


S2 


0.191 


0.170 


0.156 


0.178 





0.84 


S3 


0.255 


0.178 


0.162 


0.200 





0.85 


Rl 


0.103 


0.169 


0.158 


0.111(0.181) 


0.667 


0.48 


R2 


0.118 


0.178 


0.165 


0.120(0.194) 


0.648 


0.48 


R3 


0.136 


0.186 


0.172 


0.129(0.207) 


0.630 


0.50 



TABLE I. Central density pc, baryon rest- mass M*, ADM mass M, compactness M/R, angular 
momentum J in units of M^, and numerical results for fundamental radial oscillation period with 
/ = 2 in units of jM^I"^ of neutron stars that we pick up in this paper. Here, R denotes the 

circumference radius at the equatorial surface. For the rotating stars, the compactness measured by 
the polar radius is also listed in the bracket. All the quantities are shown in units oic = G = K = \. 



We adopted six models of neutron stars referred to as S1-S3 and R1-R3 in this numerical 
experiment. Models S1-S3 are spherical stars and R1-R3 are rigidly and rapidly rotating 
stars whose angular velocities at the equatorial surface are approximately equal to the Ke- 
plerian angular velocity (i.e., at mass-shedding limits). By exploring rotating stars at the 
mass-shedding limits, the effects of rotation on rigidly rotating neutron stars are clarified 
most efficiently. 

The maximum gravitational masses of spherical stars and rigidly rotating stars with 
n = 1 (r = 2) are ^ 0. 164/^1/2 and ^ O.lSSfsT^/^ respectively [25]. Thus, the models 
adopted here are sufficiently general relativistic in the sense that their masses are close to 
the maximum values. In Table I, characteristic quantities for models S1-S3 and R1-R3 are 
listed in the nondimensional units (in the units with c = G = K = \). 

To induce nonspherical stellar oscillations to nonrotating stars, we superimposed a ve- 
locity perturbation as 
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Su^ — —Vzu and 5uz — Vz, (29) 

where 5u^ and 5uz denote the four- velocity of cyhndrical {zu — + y^) and z components, 
y is a constant and put as V = 0.1 /zUg where zu^ denotes the coordinate radius at the 
equator: i.e., at the equatorial surface, the velocity is ~ 10% of the light speed. 

For the rotating stars, two types of the perturbations are adopted. One is a velocity 
perturbation given by 

V 

Su^ — ——w and 5uz — Vz, (30) 

with V — O.S/we. The other is a pressure perturbation in which we simply reduced the 
pressure uniformly by changing the polytropic constant from the equilibrium value to a 
smaller one. In this quasiradial oscillation is induced. Since the rotating star is 

nonspherical, gravitational waves are emitted even in this setting. 

The simulations were performed changing grid spacing and location of outer boundaries. 
As found in [12], the numerical results are sufficiently convergent if the stellar radius is 
covered by more than 30 grid points. Taking into account this fact, the grid spacing is fixed 
in the simulations of this paper as follows: For nonrotating stars, we chose it as OTe/40, and 
for rotating stars, w^j^O: Since the axial ratio of the rotating stars at the mass-shedding 
limit is ~ 0.59, the polar axis is initially covered by about 36 grid points with this setting. 

On the other hand, the location of the outer boundaries was changed varying N from 
480 to 720. We typically choose N — 720, since with this number, L is larger than a 
characteristic gravitational wave length A and, thus, the outer boundaries are located in a 
local wave zone. Since L > A, we expect that the amplitude of gravitational waves could be 
calculated within ~ 10% error. 

Table II describes the values of L and the location at which the gauge-invariant 
variables are extracted. Typically, is chosen to be ~ 0.9L. Note that varying Le from 
L/2 to 0.95L, it is found that the amplitude of gravitational waves depends very weakly on 
the location of [see Figs. 1(b) and 2(b)] for a fixed value of L. 
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Nonrotational 


L/M for N = 480, 600, 720 


Le/M for = 480, 600, 720 


X/M 


SI 


64.4, 80.5, 96.6 


53.5, 66.9, 80.3 


94.8 


S2 


57.5, 71.9, 86.3 


47.8, 59.8, 71.8 


70.5 


S3 


55.4, 69.3, 83.2 


46.1, 57.6, 69.2 


60.5 



Rotational 


L/M for = 480, 720 


Le/M for N = 480, 720 


A/M 


Rl 


63.3. 91.9 


57.9. 86.9 


(SO. 8 


R2 


58.0, 87.0 


53.1, 79.6 


73.3 


R3 


53.2, 79.8 


48.7, 73.0 


66.5 



TABLE II. Values of L and in units of M for various values of N. For comparison, we show 
a gravitational wave length for the fundamental quadrupole (Z = 2) mode derived from numerical 
results. 



B. Systematic numerical errors 

Since the outer boundaries are located in the local wave zone and, hence, the bound- 
ary conditions which are appropriate only for the distant wave zone are not precise ones, 
systematic numerical errors are generated around the outer boundaries. As a consequence, 
gravitational waveforms (gauge-invariant variables) are contaminated by numerical errors. 
To accurately extract gravitational waves, we need to eliminate the errors from raw data 
sets of the gauge-invariant variables. 

First, we summarize the behavior of the numerical errors. In Figs. 1(a) and 2(a), we 
display time evolution of raw data sets of the gauge-invariant variables extracted at several 
radii for models S2 and R2 with = 720. The gauge-invariant variables are composed 
mainly of three components; (i) a pure wave component which denotes gravitational waves, 
(ii) a constant component, and (iii) a secular drift component whose magnitude increases 
with time and is larger for the larger value of radius. To obtain clear gravitational waveforms, 
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FIG. 1. (a) Gauge-invariant variables with / = 2 in units of M for oscillation of a nonro- 
tating neutron star S2 extracted at r/M = 50.2 (dotted-dashed curve), 57.4 (dashed curve), 64.6 
(long-dashed curve), and 71.8 (solid curve), (b) gravitational waveforms after systematic numerical 
errors are subtracted. In the simulation, N = 720. 



it is necessary to eliminate the components (ii) and (iii). 

The presence of the component (ii) is in part due to the fact that the gauge-invariant 
variables are computed at finite radii. They are composed not only of gravitational waves but 
also of a quasistationary component which falls off as r~" for n> 2. In particular, rapidly 
rotating stars are of spheroidal shape and, hence, they yield quasistationary quadrupole and 
higher multipole moments which slowly vary throughout the simulation. As Fig. 2(a) shows 
(see for t — < 0), this component is smaller for the larger value of extracted radius because 
it falls off as r~" with n > 2. 

A part of the component (ii) and the component (iii) are numerical errors associated with 
an approximate treatment of the outer boundary conditions. In the following, we explain 
the origin of it in detail. 

In the wave zone, the magnitude of hij is small and, hence, it approximately obeys a 
linearized Einstein equation. In the formalism which we currently use [19], the linearized 
equation for hij is of the form 
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FIG. 2. (a) Gauge-invariant variables with / = 2 in units of M for oscillation of a rotating 
neutron star R2 extracted at r/M = 50.7 (dotted-dashed curve), 57.9 (long-dashed curve), 65.2 
(dashed curve), 72.4 (dotted curve), and 79.6 (solid curve), (b) gravitational waveforms after 
systematic numerical errors are subtracted at r/M =65.2 (dashed curve), 72.4 (dotted curve), and 
79.6 (solid curve). The simulation was performed with N = 720. 



(31) 



where Aj denotes the flat Laplacian and S^^ is composed of spatial derivatives of a = a — 1, 
(3\ p = e'^ — 1, and as 



1 



4 = (4p + 2a) - ;;V^Afi^P + 2a) 

2 „ 2 



k\l 



Here, \i denotes the covariant derivative with respect to the flat metric rjjk. 
Because of the functional form of 5*^^-, the solution of hij may be written as 

GW 2 
hij = h^j + + - -^Vij^ \kj 

where hfj^ and obey the following equations: 

e'fc = Af^k + $k + {2p + a)|fc - Fk. 



(32) 



(33) 

(34) 
(35) 



Since the inhomogeneous solution of hij, which is associated with S^j, is written by a gauge 
variable ^k, S^j does not contribute the gauge-invariant variable*. 

Gravitational waves may be extracted from G2, which is calculated by 

G2^j^f{ho-h0^W2odS. (37) 

This variable can be regarded as gravitational waves in gauge conditions in which ^fe = 
(e.g., in the harmonic gauge condition). However, in the present case, is not guaranteed 
to be vanishing and the effects of are contained in G2- Consequently, the waveforms 
may be deformed by unphysical modulations and secular drifts due to the presence of 
This illustrates that the gauge-invariant technique plays an important role for extraction of 
gravitational waves in general gauge conditions. 

General forms of outgoing-wave solutions of Eq. (34) for h'^^ have been derived by 
Teukolsky [24] and Nakamura [27]. However, numerically, such solutions can be exactly 
computed only when (A) a strict outgoing-wave boundary condition well-suited in a local 
wave zone is imposed and (B) the transverse condition is guaranteed. In the present numer- 
ical simulations, these conditions are not satisfied strictly. Therefore, unphysical modes as 
well as numerical errors contaminate numerical solutions of hf^. 

One of the candidates for the dominant unphysical modes is a solution for the equations 
Afhf^ = and /ig^ = 0. Prom the relation /ig^ = 0, /ig^ is written as Hfj{x) + Hl-{x)t 



*The gauge-invariant variables are computed from 7jj — r/y = Hij. In the linear order, H^j = 
hij + ^prjij. The linearized Einstein equation for Hij in our formulation [16,19] is 

Hij= AfHij + (4p + 2a)\ij + VikP)j + rijkP% - (Pi^j + F,-,,) - '^Viji&^fP - v'''Fk\l)- (36) 

If the Hamiltonian constraint in the linear level is satisfied, 8A/p— r/'^'F^li = 0- Then, the solution of 
Hij is hij + ^i\j + ^j\i, and consequently, the gauge-invariant variables do not contain ^j. However, if 
the Hamiltonian constraint is violated, the effects of is contained in the gauge-invariant variables. 
In the present work, we found that its effect is small, and hence we ignore it. 
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where if? and Hlj satisfy the Laplace equation. Note that H^j {n — and 1) do not satisfy 
the transverse-traceless condition in general. By performing a spherical harmonic decom- 
position of hij in the spherical polar coordinates (see Appendix A), we find the asymptotic 
behavior of the solutions of l-th moment for r » M as 

H^j r^^^, r\ r-'±\ and r'^-^. (38) 

The solutions of r^~^ and r~^~^ can be written in the form Vi\j + Vj\i where Vi denotes a 
vector. These components are eliminated from the gauge-invariant variables. However, other 
solutions cannot be written in terms of V^. Thus, four of six solutions may be contained in 
the gauge-invariant variables as unphysical modes. Properties of the unphysical solutions 
are summarized as follows: (a) they may increase linearly with time and (b) they may be 
larger for larger extracted radii because of the presence of the modes proportional to r'"'"^ 
and rK These properties agree with the numerical results shown in Figs. 1(a) and 2(a). 

Besides the global numerical errors described above, unphysical local waves of the form 
fit ± r) may be contained in hfy^. However, this is not a systematic error and, hence, it is 
not possible to eliminate systematically. Fortunately, the magnitude of such components is 
not as large as that of the systematic errors (see below) . 

The systematic non-wave components may be fitted using a function of the form C = 
Co{r) + Ci{r)t. Thus, we determine Co and Ci and then subtract C from the gauge-invariant 
variables. As mentioned above, Co(r) arises both from a nonwave component associated 
with the stationary quadrupole moment and from the unphysical modes associated with 
H^j. Ci{r)t arises from the unphysical modes associated with H^j. 

For the nonrotating stars, Co and Ci appear to be unchanged throughout the simulations. 
Thus, to determine Co(r) and Ci(r) at each radius, a least-square fitting against the gauge- 
invariant variables is carried out in the time domain. For the fitting, all the data sets with 
t — >0 are used. 

For the rotating stars, on the other hand, Cq and, in particular, Ci suddenly change 
at i — r* ~ 300M. (The reason is not very clear.) Thus, these coefficients are separately 
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determined for i — r* ^ 300M and t — r*^ 300M, carrying out the least-square fitting witfi 
two different data sets. Namely, we subtract a function of the form 

Co{r) + C,{r)t fort<t^(r), ^ ^ 

(39) 

C'oir) + C[ir)t for t>tmir), 
where tmir) is a time ~ 300M and satisfies the relation CQ{r) + Ci{r)tm{T) = Co(^) + 
Ci{r)tm{r) at each radius. 

To validate that modified waveforms depend weakly on the subtraction method, the 
following alternative method is also adopted: According to Eq. (38), Cq and Ci for are 
expressed by linear combinations of the functions r^, r^, r~^, and r~^. For a large value of 
the extraction radius Lg, we may expect that Co ~ and Ci oc at the leading order. In 
this assumption, C may be computed as 



p^r.fifW-r.j;f(r.)^,_ (40) 



where ri and r2 denote two different radii which is close to Lg. After the subtraction of the 
dominant part, it is found that a component of small magnitude associated with Cq still 
remains. To eliminate this remaining nonwave part, we simply subtract a constant from the 
resulting waveform. 

In addition to the least-square fitting method, we have also adopted this method and 
confirmed that the subtracted waveforms by this alternative method agree approximately 
with those in the least-square fitting (the wave phases agree well, and relative difference 
of the amplitude is within 10%). Thus, in the following, numerical results based on the 
subtraction using the least-square fitting are presented. 

In Figs. 1(b) and 2(b), we display improved gravitational waveforms r/i+/ sin^ 9 obtained 
after the nonwave components and the numerical errors are subtracted. It is found that the 
resulting wave amplitude and phase depend very weakly on the extracted radii for t — r* > 0. 
This confirms that the extracted components are certainly gravitational waves. 

For the rotating star, the amplitudes of the gravitational waveforms extracted at different 
radii are in slight disagreement each other. Figure 2(b) shows that the magnitude of the 
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difference is ~ 10~^. This indicates that even in the improved waveforms, numerical errors 
of magnitude ~ 10~^ still remain. The origin is likely to be a nonsystematic error such as 
spurious wave components. 

Besides quasiperiodic waves, a spike is visible at the beginning of the simulation (at 
t ~ lOM) in Figs. 1 and 2. We do not understand the origin of it exactly. The following is 
inference for the possible origin: At t=0, we rather crudely add a weakly nonlinear velocity 
perturbation to equilibrium states. Because of the weak nonlinearity, impulsive gravitational 
waves may be excited besides eigen oscillation modes of the equilibrium stars. Such impulsive 
gravitational waves seem to propagate at t — r ~ 0. 

Before closing this subsection, we note the following: The magnitude of the numerical 
error associated with H^^ is not very large in three-dimensional simulations for binary neu- 
tron star merger [19] and oscillation of neutron stars [17], although it might be contained. 
We suspect that the excitation of such unphysical modes may be associated with the in- 
terpolation used in the Cartoon method. (In this paper, we simply adopt a second-order 
interpolation.) This suggests that there would be still a room to improve the interpolation 
scheme in our numerical implementation. 

C. Gravitational waveforms 

1. Nonrotating stars 

Improved gravitational waveforms with / = 2 from axisymmetrically oscillating and non- 
rotating neutron stars for models S2 and S3 are displayed in Fig. 3. For both models, 
the simulations were performed with = 480 (dotted-dashed curves), 600 (long-dashed 
curves), and 720 (solid curves). Gravitational waveforms evaluated by the quadrupole for- 
mula (dashed curves) are displayed together to illustrate its validity. The quadrupole formula 
was used in all the simulations and it is found that the gravitational waveforms depend very 
weakly on the value of N. Here, the numerical results for N = 720 are plotted. 
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FIG. 3. Gravitational waves with I = 2 from nonrotating neutron stars with axisymmetric 
oscillations (a) for models S2 (left) and (b) S3 (right). For both models, we show the results 
with = 480 (dotted-dashed curve), 600 (long-dashed curve), and 720 (solid curve). The dashed 
curves denote the corresponding gravitational waveforms by the quadrupole formula. Here, the 
gravitational waveforms are extracted at ^ 0.9L. Note that the behavior of the raw data is 
very similar to that in Fig. 1, i.e., with the increase of Lg (with the increase of for a fixed grid 
spacing), the amplitude of the drift is larger. 

Figure 3 indicates that one oscillation mode is dominantly excited. As a result, the 
waveforms are well-approximated by a simple sine curve. Indeed, the Fourier spectra of 
^^quad^ sin^ 6 posscss oue sharp peak, and the oscillation periods are ~ 0.84, 0.84, and 
0.85 in units of 2tt ^R-^/M for models SI, S2, and S3, respectively ^ Thus, irrespective 
of compactness of the neutron stars, the oscillation period is ~ 0.85 x 2ti^JbFJm. This 
implies that the oscillation is associated with the fundamental quadrupole mode, since the 



^In [15], a three-dimensional simulation for an oscillating and nonrotating star which is almost the 
same as model SI was performed. They computed approximate gravitational waveforms in a near 
zone (A/L 0.12) and derived the period for the fundamental quadrupole mode as 0.83-^"^/ K-^/M. 
This value agrees with ours in ~ 1% error. 
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coefficient (~ 0.85) depends very weakly on tlie compactness of tlie neutron stars for a given 
equation of state [26] . 

From Eq. (14), the luminosity of gravitational radiation as a function of time is com- 
puted. Since the luminosity also varies in a periodic manner, we define an averaged energy 
flux according to 



where to is a constant. For models S1-S3, the averaged luminosity is ~ 6 x 
10~® (in units of G~^c^). Therefore, the energy dissipated in one oscillation period is much 
smaller than the total mass energy of the system and the radiation reaction timescale is 
much longer than the oscillation period. 

Figures 1 and 3 indicate that (i) the wave length is independent of L and Lg, (ii) the 
amplitude of gravitational waves is overestimated for the smaller values of L and Lg, and 
(iii) the amplitude for model S3 approximately converges to an asymptotic value for N > 
600 within ~ 10% error. These facts suggest that for Lg ^ A (see Table II), convergent 
gravitational waveforms within 10% error can be computed. On the other hand, if the outer 
boundaries are located in a near zone with L < X, the amplitude of gravitational waves 
is overestimated: For L pa 2A/3, it is overestimated by a factor of ~ 20%. This result is 
consistent with that in our previous study for gravitational waves from binary neutron stars 
in quasiequilibrium circular orbits [28]. 

The quadrupole formula yields well-approximated gravitational waveforms besides a sys- 
tematic underestimation of the wave amplitude. For models S2 and S3, the asymptotic 
amplitudes of A2/M are about 0.010 and 0.007, respectively. On the other hand, according 
to the quadrupole formula, they are about 0.008 and 0.0055. Thus, the underestimation 
factor is ~ 20%. If it is due to the first post Newtonian correction, the factor should be 
proportional to the compactness of neutron stars M/R(— GM/Rc^) or v'^{— jc?) where 
V denotes typical magnitude of the oscillation velocity. To determine which the dominant 
component is, we performed a simulation reducing the magnitude of the velocity perturba- 




(41) 
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tion initially given (i.e., reducing the magnitude of V) and found that the wave amplitude is 
underestimated by ~ 20% irrespective of the magnitude of V. Therefore, we conclude that 
the underestimation factor is proportional to the magnitude of M/ R in the present case. 

Although the wave amplitude is underestimated, the wave phase is computed accurately 
by the quadrupole formula. The most important element in detection of gravitational waves 
using matched filtering techniques is to a priori know the phase of gravitational waves. From 
this point of view, the quadrupole formula is a useful tool for computation of gravitational 
waveforms. 

2. Rotating stars 

In Fig. 4, gravitational waveforms with I — 2 from oscillating and rapidly rotating neu- 
tron stars for models R2 and R3 are displayed. For these simulations, velocity perturbations 
are added initially without changing other quantities. The numerical results are shown for 
N = 480 (long-dashed curve) and 720 (solid curve). For N = 480 and 720, the gauge- 
invariant variables are extracted at ~ 2 A/3 and A, respectively. The waveforms computed 
by the quadrupole formula (dashed curves) are displayed together. 

In contrast to those from oscillating and nonrotating neutron stars, the gravitational 
waveforms are not of simple sine curve. There are two main reasons. One is the following: 
In the nonrotating case, the restoring forces against a compression for zu and z directions 
are of identical magnitude and, hence, the oscillation periods of two directions agree. For 
the rotating stars, on the other hand, the two restoring forces are not of identical magnitude 
and, therefore, the oscillation periods of two directions do not agree. Due to this fact, the 
waveforms are composed of more than two oscillation modes. 

The other reason is that gravitational waves are emitted due to a quasiradial motion in 
the case of rotating neutron stars in contrast to nonrotating neutron stars. In particular, 
we here choose rapidly rotating neutron stars and, therefore, the wave amplitude can be as 
larger as than for the quadrupole oscillations. 
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FIG. 4. Gravitational waves with / = 2 from rotating neutron stars R2 (left) and R3 (right) 
with N = 480 (long-dashed curve) and 720 (solid curve). The dashed curves denote the results by 
the quadrupole formula. 



To clarify what oscillation modes are relevant, the Fourier spectrum F{f) for models 
R2 and R3 are displayed in Fig. 5. This shows that there are two characteristic peaks in 
the spectrum. The oscillation periods defined by l//peak where /peak is the frequency of the 
peaks in the Fourier spectrum are listed in Tables I and III. 

For models R1-R3, the oscillation period of the larger frequency is 

Pose ~ 0.5 x27r^ — , (42) 

where R is the equatorial circumferential radius. The coefficient (~ 0.5) depends very 
weakly on the compactness of the neutron star. This indicates that the oscillation mode is 
the fundamental quadrupole mode. 

The frequency of the other peak is smaller than that of the fundamental quadrupole 
mode. This peak is likely to be associated with the fundamental quasiradial oscillation 
mode {pi mode). To confirm this prediction, we performed the Fourier analysis to the 
time sequence of the central density and found that the characteristic oscillation period 
indeed coincides with that of the second peak (see Table III). Furthermore, it agrees with 
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FIG. 5. The Fourier spectrum of gravitational waveforms for models R2 (left) and R3 (right). 

The spectrum is normalized by the maximum value, and the frequency / is shown in units of M~^. 

The peaks of the smaller and larger frequencies indicate the quasiradial and quadrupole modes, 

respectively. 

the characteristic frequency for quasiradially oscillating neutron stars presented in [12]. All 
these facts confirm our prediction. 

The Fourier spectra indicate that gravitational waves from axisymmetric global oscil- 
lations of rigidly rotating stars are composed of two dominant modes. To confirm this 
conjecture, we performed simulations initiated from another initial condition in which the 
pressure is uniformly depleted but the velocity field is unperturbed. To uniformly deplete 
the pressure, K is decreased by 20% at t — 0. In Fig. 6, we display (a) the gravitational 
waveform and (b) the Fourier spectrum for model Rl for this simulation. Figure 6 (b) shows 
that two modes, the fundamental quadrupole and quasiradial ones, are dominant again. 
This result justifies our conjecture. In contrast to the cases in which the nonspherical ve- 
locity perturbation is added, the mode with lower frequency, i.e., the quasiradial mode, is 
dominant. This is because the matter motion is almost quasiradial. 

In the collapse of massive rotating stellar core, a protoneutron star will be formed. If 
the progenitor star is not rapidly rotating and its degree of differential rotation is not high. 
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the protoneutron star relaxes to a nearly quasistationary state soon after the collapse (e.g., 
[1]). At the formation of a rotating protoneutron star in a nearly quasistationary state, 
nonspherical oscillations are excited by the quasiradial infall. Because of the nonspherical 
nature, gravitational waves are emitted [1]. As illustrated above, in such oscillating neutron 
stars in a nearly quasistationary state, two dominant modes (quasiradial and quadrupole 
modes) may be excited. If the collapse is quasiradial, the quasiradial mode will be the main 
component. On the other hand, if the nonspherical quadrupole oscillation of high amplitude 
is excited, the quadrupole mode will be dominant. 




(a) (b) 

FIG. 6. (a) Gravitational waves with / = 2 from rotating neutron stars Rl in a quasiradial 

oscillation with = 480 (solid curve). The dashed curves denote the results by the quadrupole 

formula, (b) The Fourier spectrum of gravitational waves. The spectrum is normalized by the 

maximum value, and the frequency / is shown in units of M~^. The peaks of the smaller and 

larger frequencies indicate the quasiradial and quadrupole modes, respectively. In this case, the 

quasiradial mode is dominant (compare with Fig. 5). 



As Fig. 4 shows, the wave amplitude decreases with time in the early phase (t — r^, ^ 
500M), and then relaxes to a small value. On the other hand, the amplitude does not de- 
crease in Fig. 6. This indicates that the amplitude of the quadrupole oscillation decreases 
with time to a small value, while the quasiradial oscillation does not. As mentioned above, 
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the oscillation frequencies for z and w directions are not identical for the quadrupole oscilla- 
tion of rotating stars. As a result, shocks may be formed at a collision of compression waves 
in the two different oscillation directions. Repeating this process, the quadrupole oscillation 
may be gradually damped. On the other hand, there is no process that damps quasira- 
dial mode quickly. Thus, it is reasonable to expect that the quasiradial mode eventually 
dominates in the oscillating and rotating stars. 

Using Eq. (14), the luminosity is computed. It is found that until t ~ lOOOM, the 
total radiated energy is computed as ~ 5 x 10~^M, which is much smaller than the total 
mass energy of the system. Therefore, the damping timescale of the wave amplitude due to 
gravitational waves is much longer than the oscillation and rotation periods. 

As in the nonrotating case, approximate gravitational waveforms were computed using 
the quadrupole formula. Figures 4 and 6 indicate that the waveforms agree with those 
computed by the gauge-invariant method, besides a systematic underestimation of the am- 
phtude. The underestimation factor is of order Mj R as in the nonrotating cases. In the case 
of rotating stars, a modulation of the wave amplitude is outstanding. Using the quadrupole 
formula, however, such modulation can be well computed. As mentioned in IV C 1, the most 
important element in detection of gravitational waves using matched filtering techniques is 
to a priori know the phase and modulation of gravitational waves. The results here indi- 
cates that they are computed well in the quadrupole formula. Therefore, for computation of 
gravitational waves in rotating stellar core collapse to a protoneutron star, the quadrupole 
formula may be a useful tool. 

V. SUMMARY AND DISCUSSION 

We have studied gravitational waves from axisymmetrically oscillating neutron stars 
adopting the gauge-invariant wave extraction method as well as the quadrupole formula. 
It is found that several types of the nonwave components such as the stationary parts of 
metric and numerical errors are contained in the gauge-invariant variables. The numerical 



29 





Posc/(27ri?^/^M"^/^)(for central density) 


-fosc/(2vri?'^/^-M ^/^)(for gravitational waves) 


Ref. [12] 


Rl 


0.66 


0.66 


0.66 


R2 


0.69 


0.70 




R3 


0.76 


0.76 


0.76 



TABLE III. Posc/(27ri?3/2M-V2) for the fundamental quasi-radial mode calculated from the 
evolution of the central density and gravitational waveforms. The last column shows the results 
obtained from simulations of the quasiradial oscillation. 



errors are generated due to an approximate treatment for the outer boundary conditions. 
We illustrate a method to subtract the dominant components of the numerical errors and 
demonstrate it possible to extract gravitational waves even from such noisy data sets with 
a residual of magnitude ~ 10~^. 

The gravitational waveforms computed in the quadrupole formula agree well with those 
obtained from the gauge-invariant technique besides a systematic underestimation of the 
amplitude by ~ 20%. An important point is that the evolution of the wave phase and the 
modulation of the amplitude are computed with a good accuracy. This indicates that for a 
study of gravitational waveforms from rotating stellar core collapse to a protoneutron star, 
the quadrupole formula will be a useful tool in fully relativistic simulations. It should be also 
addressed that the result in this paper supports the treatment in [1] in which gravitational 
waveforms are computed using a quadrupole formula in approximate general relativistic 
simulations. 

The gauge-invariant variables are extracted for various values of extraction radii. It is 
found that to extract gravitational waves within ~ 10% error, the extraction radius has to 
be larger than ~ 90% of the gravitational wave length. If the outer boundaries are located 
in the near zone with L < A, the amplitude of gravitational waves is overestimated: For 
L fa 2 A/ 3, it is overestimated by ~ 20%. For L < 2 A/ 3, the factor of the overestimation is 
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even larger. 

In the present work, the amplitude of gravitational waves in a local wave zone is much 
larger than that of systematic numerical errors. This fact enables to subtract them from 
the gauge-invariant variables accurately. If the magnitude of the errors is much larger than 
that of the amplitude of gravitational waves, however, it would not be possible to carry out 
an accurate subtraction. For example, in rotating stellar core collapse, the amplitude of 
gravitational waves in the local wave zone at r ~ A is at most ~ 10~^ according to gravita- 
tional waveforms calculated by a quadrupole formula [1]. To extract gravitational waves of 
such small amplitude, it is necessary to reduce the magnitude of the numerical errors. To 
achieve that, we need to impose more accurate outer boundary conditions. Developing such 
conditions is crucial in computing gravitational waves of small amplitude of O(10~^) from 
raw data sets of metric. 

Another possible method for computing accurate gravitational waves of small amplitude 
is to adopt a quadrupole formula taking into account higher-order post Newtonian terms. 
As indicated in this paper, the simple quadrupole formula underestimates the amplitude of 
gravitational waves by 0{M/R). In rotating stellar core collapse, the error in the amph- 
tude will be ~ 10%. To compute the amplitude within ~ 1% error, we should take into 
account higher general relativistic corrections. In quadrupole formulas with the higher post 
Newtonian corrections (as derived in [29]), it may be possible to obtain gravitational wave- 
forms within 1% error. Such formulas will be useful to extract gravitational waves of small 
amphtude from rotating stellar core collapses and from oscillating neutron stars. 

In addition to the study for gravitational wave extraction, oscillation modes of rotating 
neutron stars are analyzed. It is found that two modes (the fundamental quadrupole and 
quasiradial modes) are dominantly excited due to the global oscillation. The frequency of the 
quadrupole mode is proportional to and is higher than that of the quasiradial one 

for the typical values of mass and radius of neutron stars. It is shown that the amplitude of 
the quadrupole mode decreases with time due to an incoherent nature of the oscillation, but 
that of the quasiradial mode is not damped quickly, hence being the dominant mode after 
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several dynamical timescales. We expect that in rotating stellar core collapse to a protoneu- 
tron star in a nearly quasistationary state, these two modes may be the main components 
in the burst phase of gravitational waves. The quadrupole mode will be damped within 
a few dynamical timescales and subsequently the quasiradial mode will be the dominant 
component to be longterm quasiperiodic waves. 
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APPENDIX A: SOLUTIONS FOR THE TENSOR LAPLACE EQUATION 

Here, we describe solutions for the tensor Laplace equations in spherical polar coordinates 

as 



Afhab = 0. 



(Al) 



hab is expanded by tensor harmonic functions as 



hab = J2 



( 



+ 



AiYi^ rBiYio,e 
* r^KiYio + GiWio) 

r^sm^e{KiYio-GiWi 

\ 



* * 



10 



rCideYio sin 9 

* -r^DiWio sin 9 

* * 



(A2) 



where Ai, Bi, Ci, Di, Gi, and Ki are functions of r. With the above expansion, the compo- 
nents of the Laplace equation are written as 
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Afh 
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1/0 = 0, 



deYi^ = 0, 
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+ -^i —^l H o 



9 * 9 



Afh, 



I '- ' 



Xi-2 



Im 
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sin 6* 
sin eWio = 0, 
2Ai 



= 0, 



Ki + -Ai - -^Bi ]Yio+[ F!' + -F: 



-Fi^-BAWi 



'10 



0, 



sin^ 9 



J2(-H^Yio-HlWi 



10 



0, 



(A3) 



where A; = /(/ + 1) and ' denotes d/dr. 

Setting Ai — Air'^, Bi — Bir'^, Ki — Kir'^, and Fi — Fir'^ where Ai ~ Fi are constants, 
simultaneous equations for the even-parity modes are derived as 



X - 4 4Ai 4 

2 x-4 -2 2Ai-4 

2 -2Ai x-2 Q 

2 x+2 



Ai 
Bi 

K, 



yF,j 



0, 



(A4) 



where x = n{n + 1) — A/. For the existence of the solutions, the determinant of the matrix 
should be zero. Then, an algebraic equation for n is derived, and the solutions are n — l±2, 
I, and —I — 3. The relations among Ai ~ Fi are easily obtained for each value of n. 

It is also easy to check that the solutions with n = I — 2 and — / — 3 are written in the form 
Vi\j + Vj\i using a vector of even-parity. 

Prom the same procedure, the solutions for the odd-parity mode are written as Q = Qr" 
and Di = Dir^ where n = / ±1, — / and —1^2, and Ci and Di are constants. In this case, the 
solutions with n — I — 1 and —1 — 2 are written as Vi\j + Vj\i using a vector of odd-parity. 
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